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Introduction 


This  final  report  documents  a  Phase  I  SBIR  project  to  implement  and 
conduct  initial  testing  of  an  algorithm  to  perform  multisensor  tracking  of  an 
aircraft  based  on  the  fusion  of  position  measurements,  angle  measurements, 
and  optical  images  of  the  aircraft.  The  work  was  performed  by 
Wagner  Associates  under  contract  to  Rome  Laboratory,  with  consulting 
support  from  Dr.  Michael  Miller  of  the  Department  of  Electrical  Engineering 
at  Washington  University  in  St.  Louis  (WUSL). 

The  principal  focus  of  the  Phase  I  study  was  to  adapt  to  a  serial 
computing  workstation  environment  (Tatung  Super  COMPstation  20) 
algorithmic  procedures  for  executing  automatic  target  detection,  tracking,  and 
recognition  (ADTR)  which  were  developed  at  WUSL  and  implemented 
across  a  network  of  Silicon  Graphics  workstations  and  a  massively  parallel 
DECmpp  12000/Sx  (8192  processors). 

We  proceed  to  describe  the  general  model  framework  of  our  analysis. 


Sensors 

The  operational  setting  we  assume  in  Phase  I  is  that  of  a  single  aircraft 
of  known  type  being  tracked  in  combination  by  a  low  resolution  radar  which 
provides  point  estimates  of  the  target  location  and  a  high  resolution  optical 
sensor  which  provides  (i)  measurements  of  the  direction  to  the  target  and 
(ii)  gray  scale  pixel  images  of  the  target  body. 

We  model  the  low  resolution  sensor  as  providing  relatively  inaccurate 
measurements  of  target  location  (cJr  =  200  meters,  =  0.5°,  0^  =  2°).  These 
measurement  errors  are  intended  to  be  representative  of  military  air  traffic 
control  radars  and  shipboard  surveillance  radars. 


The  high  resolution  sensor  provides  more  accurate  target  azimuth  and 
elevation  angle  measurements  (0^2=00=0.!°)  than  the  low  resolution 
sensor  but  provides  no  direct  range  measurements  In  addition,  the  high 
resolution  sensor  is  assumed  to  generate  focal  plane  pixel  images  of  the  target. 
A  library  of  reference  image  templates  has  been  provided  to  us  by  WUSL  for 
an  optical  system.  Each  template  in  the  library  is  a  64  x  64  pixel  array  with 
each  pixel  representing  a  signal  intensity  level  discretized  on  a 
monochromatic  (gray  scale)  of  256  possible  levels. 


Target  Tracking 

In  the  present  context,  a  target  track  is  the  joint  specification  of  the 
position,  velocity,  and  orientation  history  of  the  target  over  a  time  window 
(possibly  a  moving  window)  and  of  the  target  aircraft  type. 

The  ADTR  algorithm  automates  the  track  estimation  function  by  first 
applying  classical  Kalman  methods  to  process  the  low  resolution  sensor  data. 
The  Kalman  filter  tracking  solution  is  used  to  initialize  the  high  resolution 
tracking. 

The  target  image  data  provided  by  the  high  resolution  sensor  is 
processed  in  a  quite  different  manner.  In  this  case,  a  Monte  Carlo  scheme  is 
used  for  random  sampling  from  the  Bayes  posterior  distribution  on  the  joint 
position/velocity/orientation  state  of  the  target.  The  reported  mean  target 
state  estimate  and  the  associated  error  covariance  are  based  on  the  sample 
statistics  of  this  posterior  distribution. 

A  principal  objective  of  our  Phase  I  study  is  to  quantify  the  reduction  in 
tracking  error  that  derives  from  the  fusion  of  the  high  resolution  sensor  data 
with  the  low  resolution  sensor  data.  There  are  two  important  effects  at  work 
here.  First,  the  more  accurate  angle  measurements  provided  by  the  high 
resolution  sensor  accelerate  track  solution  convergence  and  reduce  tracking 
error.  Second,  estimates  of  the  target  orientation  angles  (in  particular,  pitch 
and  yaw)  derived  from  the  processing  of  image  data  can  be  used  to  draw 
inference  about  the  instantaneous  direction  of  target  flight.  Thus,  there  is  the 
possibility  of  the  early  recognition  of  directional  changes  by  the  target.  If  only 
target  position  can  be  observed,  the  occurrence  of  a  maneuver  must  be 
inferred  through  the  measurement  of  subsequent  target  positions.  The 
recognition  of  a  target  maneuver  can  then  significantly  lag  its  occurrence. 

As  noted,  in  Phase  I  we  address  only  the  single  target  tracking  problem. 
In  Phase  n  we  plan  to  extend  the  ADTR  algorithm  to  the  multitarget  tracking 
problem.  We  then  expect  to  observe  further  benefits  of  multisensor  fusion. 
Multiple  targets  are  inherently  easier  to  separate  in  problem  space  when 
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aircraft  orientation  and  aircraft  classification  are  included  as  state  variables. 
Improved  separation  translates  into  a  reduced  incidence  of  report-to-track 
association  errors  and  therefore  into  increased  tracking  accuracy. 


Core  Algorithm  Design 

In  Chapter  I  we  document  the  methods  that  underlie  the  ADTR 
algorithm.  The  basic  approach  is  to  use  Monte  Carlo  simulation  methods  in 
conjunction  with  a  random  sampling  procedure  to  construct  sample  target 
tracks  which  incorporate  the  fused  location  and  image  data  provided  by  the 
sensors. 

The  sampling  procedure  we  employ  is  based  on  the  jump-diffusion 
model  developed  by  Grenander  and  Miller  (see  reference  [1])  and  adapted  in 
references  [2]  and  [3]  to  the  aircraft  ADTR  problem.  The  key  idea  is  to  allow  a 
stochastic  differential  equation  (SDE)  and  a  Markov  jump  process  to  act  in 
concert  on  the  Monte  Carlo  samples  of  the  target  track.  The  role  of  die  SDE  is 
to  act  on  the  continuous  target  state  variables  (position,  velocity,  and 
orientation)  and  smoothly  deform  the  Monte  Carlo  sample  paths  into  sample 
tracks  drawn  from  the  Bayes  posterior.  The  role  of  the  jump  process  is  to  act 
on  the  single  discrete  target  state  variable  (aircraft  type)  to  allow  new  image 
data  to  alter  the  aircraft  type  assigned  to  an  individual  Monte  Carlo  track. 
(Note:  Because  we  assume  in  Phase  I  that  target  type  is  known,  we  have 
deferred  to  Phase  II  the  implementation  of  the  jump  feature  of  the 
jump-diffusion  model.) 

The  SDE  model  (see  equation  (T3))  effectively  executes  a  randomized 
gradient  descent  search  such  that  the  term  corresponding  to  the  negative 
gradient  of  the  Gibbs  potential  function  drives  the  Monte  Carlo  paths  in  the 
direction  of  decreasing  potential  and  hence  increasing  likelihood.  The  white 
noise  term  in  the  SDE  adds  the  correct  degree  of  randomness  to  the  sampling 
process  to  generate  the  Bayes  posterior  distribution. 

The  underlying  theory  guarantees  that  the  combined  action  of  smooth 
target  position,  velocity,  and  orientation  deformation  by  the  SDE  and  the 
discrete  jump  process  for  assigning  aircraft  type  will  result  in  Monte  Carlo 
tracks  that  closely  approximate  tracks  drawn  from  the  Bayes  posterior.  An 
essential  computational  advantage  of  this  sampling  procedure  is  that  it 
accomplishes  the  sampling  without  requiring  the  computationally 
prohibitive  task  of  constructing  the  posterior  target  track  density  function. 
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Numerical  Studies 

In  Chapter  n  we  summarize  the  numerical  results  we  obtained  based 
on  our  Phase  I  implementation  of  the  ADTR  algorithm.  Our  Phase  I  testing 
is  focused  exclusively  on  parametric  studies  to  investigate  the  robustness  of 
the  tracking  function  of  the  algorithm  as  well  as  its  sensitivity  to  certain  key 
problem  parameters  related  to  the  characteristics  of  the  high  resolution  sensor 
and  the  characteristics  of  the  templates  in  the  image  library.  The  following 
table  identifies  which  problem  parameters  are  varied  and  the  assumed  base 
case  value  of  each. 


Parameter 

High  resolution  sensor 
angular  measurement  error 

High  resolution  sensor  data  rate 

Pixel  resolution 

Gray  scale  intensity  discretization 
Pixel  noise 


Base  Case  Value 

0.1° 

2Hz 

64  x64  pixels 
256  levels 
12  gray  scale  units 


The  ground  truth  target  track  assumes  an  X29  aircraft  initially  in  level 
flight  at  10,000  ft.  The  aircraft  then  executes  a  banked  (roll  angle  =  -45°)  and 
pitched  (pitch  angle  =  20°)  turn  of  180°  with  a  yaw  rate  of  5°/sec.  The 
maneuver  requires  36  seconds  to  complete. 

We  calculate  and  plot  tracking  error  (measured  as  the  distance  between 
the  estimated  target  position  and  the  actual  target  position)  versus  time  imder 
three  different  assumptions  about  the  available  sensor  data: 


(LR):  Low  resolution  position  measurements  only 

(HRA):  Low  resolution  position  measurements  +  high 

resolution  angle  measurements 

(HRI):  Low  resolution  position  measurements  +  high 

resolution  angle  measurements  +  high  resolution 
pixel  images. 
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These  tracking  accuracy  comparisons  are  shown  in  Figures  II-2 
through  n-11. 


We  summarize  our  base  case  numerical  results  as  follows: 


•  Case  I: 

(i)  target  not  maneuvering 

(ii)  fusion  of  low  resolution  (radar)  and  high  resolution 
(optical)  point  tracking  data  (HRA  processing) 

(iii)  Kalman  filter  state  estimation  methods. 

Then  tracking  errors  vary  from  200  meters  to  500  meters.  This 
compares  with  tracking  errors  of  as  much  as  3,000  meters  under 
low  resolution  tracking  (LR  processing). 


•  Case  II: 

(i)  target  maneuvering  (2g  acceleration  for  36  secs.) 

(ii)  fusion  of  low  resolution  point  tracking  data,  high 
resolution  point  tracking  data,  and  imaging  data  (HRI 
processing) 

(iii)  jump-diffusion  method  of  state  estimation. 

Then  tracking  errors  vary  from  200  meters  to  600  meters.  This 
compares  with  tracking  errors  of  as  much  as  1,400  meters  under 
high  resolution  point  tracking  (HRA  processing). 


Phase  I  Study  Conclusions 

Our  principal  study  conclusions  are: 

•  We  successfully  demonstrated  that  the  computations  required  to 
carry  out  the  data  fusion  algorithm  operations  using  our 
estimation  methods,  previously  executed  only  within  the 
context  of  a  massively  parallel  computer  architecture,  can  be 
performed  on  a  serial  processing  workstation  (the  equivalent  of 
a  Sun  Microsystems  SPARC  20). 


•  Our  study  shows  that  there  are  significant  benefits  to  be  gained 
in  terms  of  increased  tracking  accuracy  through  the  fusion  of 
multisensor  point  tracking  and  imaging  data.  We  have 
demonstrated  this  specifically  for  the  combination  of  a  low 
resolution  radar  sensor  and  a  high  resolution  optical  imaging 
sensor.  However,  the  methods  we  have  developed  to  fuse 
multisensor  data  are  quite  general  and  can  be  extended  to  other 
sensing  technologies. 


Phase  II  Plan 

In  Phase  II  we  will  extend  our  analysis  to  the  multiple  target  case.  In 
that  context  the  ADTR  function  will  be  extended  to  include  data  association, 
i.e.,  the  formation,  scoring,  and  pruning  of  hypotheses  regarding  which 
subsets  of  measurements  to  fuse  into  target  tracks.  The  jump-diffusion 
method  applied  in  Phase  I  to  target  type  determination  will  be  adapted  to 
perform  the  critical  data  association  operations. 

Our  planned  Phase  n  technical  objectives  include: 

(1)  Extend  the  application  of  the  algorithm  to  the  case  of  multiple 
air  targets. 

(2)  Expand  the  image  template  library  to  include  multiple  aircraft 
types  and  implement  the  program  logic  to  perform  the 
automatic  target  recognition  function  when  the  target  type  is 
unknown. 

(3)  Allow  for  the  aircraft  images  (both  those  stored  in  the  template 
library  and  those  observed  by  the  optical  sensor)  to  depend  on 
visibility  conditions  and  the  level  of  ambient  light. 

(4)  Incorporate  program  logic  to  deal  with  unrecognizable  targets.  A 
given  target  may  be  imrecognizable  for  a  number  of  reasons: 

•  The  target  is  an  aircraft  of  a  type  not  included  in  the  image 
template  library 

•  Multiple  targets  are  closely  spaced  causing  some  degree  of 
image  occlusion 

•  Camouflage  is  being  employed  in  an  attempt  to  defeat  a 
tracking  system  based  on  image  recognition  methods. 

(5)  Extend  the  application  of  the  algorithm  to  microwave  and 
millimeter  wave  radar  imaging  technologies  including  SAR  and 
ISAR. 
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Chapter  I 


Automatic  Detection.  Tracking  and  Recognition  (ADTR) 

Algorithm  Description 


In  this  chapter  we  document  the  algorithm  to  perform  the  automatic 
detection,  tracking,  and  recognition  (ADTR)  function  which  we  have 
implemented  as  part  of  our  Phase  I  study.  Mathematical  details  are  included 
to  the  extent  that  they  do  not  already  appear  in  published  papers  and  reports 
(see,  for  example,  references  [3],  [4],  and  [5])  and  to  the  extent  that  such  details 
contribute  to  an  understanding  of  how  the  underlying  theory  has  been 
applied  in  the  specific  context  of  our  Phase  I  study. 

As  noted  in  the  Introduction,  the  focus  of  our  Phase  I  study  has  been 
the  implementation  and  initial  testing  of  an  algorithm  to  execute  the 
function  of  the  automatic  recognition  of  an  air  target.  We  begin  in  Section  I-l 
with  a  formal  statement  of  the  ADTR  problem  addressed  in  Phase  I.  In 
Section  1-2  we  describe  the  random  sampling  technique  based  on  the 
jump-diffusion  process  model  which  underlies  our  ADTR  algorithm.  The 
specific  ADTR  algorithm  procedures  are  detailed  in  Section  1-3. 


I-l  Phase  I  ADTR  Problem  Description 

In  defining  the  automatic  target  recognition  problem,  we  distinguish 
between  those  components  of  the  target  state  vector  which  vary  continuously 
and  those  components  which  vary  discretely.  In  the  context  of  our  Phase  I 
study,  the  continuous  state  components  are  taken  to  be  the  aircraft  position 
vector  p  =  ip^,P2fP3)r  the  aircraft  velocity  vector  v-(vi,V2,V3),  and  the  aircraft 
orientation  (Euler)  angle  vector  (l)  =  (^i,(|)2/4>3)/  specified  in  an  inertial 
coordinate  frame.  Here,  <{>1  =  roll  angle,  <|)2  =  pitch  angle,  and  <])3  =  yaw  angle. 
The  single  discrete  state  variable  is  aircraft  type,  denoted  by  aeA,  where  A  is 
a  known  alphabet  of  possible  aircraft  types. 
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We  postulate  two  types  of  sensors  that  provide  information  about  the 
target  state  vector.  First,  there  is  a  low  resolution  tracking  sensor  that 
provides  measurements  of  target  position.  We  model  the  low  resolution 
sensor  as  generic  in  that  the  particular  sensing  technology  it  employs  is 
unimportant.  The  only  essential  requirement  is  that  each  sensor 
measurement  provide  an  estimate  of  target  azimuth  angle,  elevation  angle; 
and  range.  Measurement  errors  in  the  three  coordinates  are  assumed  to  be 
statistically  independent.  Possibilities  for  the  low  resolution  sensor  are 
(i)  a  narrow  band  passive  cross-array  of  isotropic  sensor  elements 
(see  reference  [3])  and  (ii)  a  tracking  radar. 

Second,  there  is  a  high  resolution  imaging  sensor  that  provides  pixel 
images  of  the  target  such  that  the  signal  in  each  pixel  is  a  noisy  gray  scale 
intensity  level.  The  case  we  treat  in  Phase  I  is  that  of  an  optical  sensor  system 
that  produces  2-D  far-field  orthographic  projections  of  the  target  image  on  the 
sensor  focal  plane.  Each  projected  image  is  represented  by  a  64  x  64  array  of 
gray  scale  intensity  levels  (discretized  as  integer  values  on  the  interval  [0,255]). 
Thus,  a  high  resolution  measurement  is  an  element  of  the  space  [0,255]^’'^. 

In  our  Phase  I  problem  formulation,  we  assume  that  the  low  and  high 
resolution  sensors  are  collocated  at  a  single  ground  tracking  station.  The 
extension  to  sensors  at  multiple  geographic  locations  is  deferred  to  Phase  n. 

The  ADTR  problem,  as  we  have  posed  it  in  Phase  I,  then  is  to  develop 
an  algorithmic  procedure  to  detect  a  single  air  target,  identify  its  type  (i.e., 
classify  the  target),  and  generate  estimates  for  its  position,  velocity,  and 
orientation.  This  combination  of  detection,  tracking,  and  recognition 
functions  is  to  be  performed  on  the  basis  of  a  combination  of  low  resolution 
and  high  resolution  sensor  data,  together  with  prior  knowledge  about  the 
candidate  aircraft  types  and  flight  characteristics. 


1-2  The  Jump-Diffusion  Method 

The  approach  we  take  to  the  ADTR  problem  is  based  on  the 
jump-diffusion  method  for  sampling  from  the  Bayes’  posterior  density  for 
the  target  track  given  a  combination  of  prior  information  about  the  target  and 
accumulated  sensor  data.  In  the  current  context,  a  target  track  X(f)  is  a 
specification  of  the  target  path  through  position,  velocity,  and  orientation 
space  at  a  sequence  of  discrete  times  f*  over  the  time  interval  [0,f]  as  well  as  a 
specification  of  the  target  typ>e.  We  take  the  to  be  the  sensor  measurement 
times. 
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1-2.1  Theory.  As  we  describe  below,  a  jump-diffusion  process  has  the 
defining  properties  that  it  (i)  executes  jumps  in  the  discrete  state  variables  at 
exponentially  distributed  times  and  (ii)  follows  a  stochastic  differential 
equation  diffusion  model  in  the  continuous  state  variables  between  the 
jumps.  The  key  theoretical  result  on  which  our  state  estimation  methods  are 
based  is  an  ergodic  theorem  proven  in  reference  [2]  which  insures  that  the 
jump-diffusion  process  we  construct  converges  in  distribution  to  the  Bayes' 
posterior  distribution  for  the  target  track  conditioned  on  all  past  information 
G?oth  prior  information  and  sensor  measurements).  This  result  provides  the 
basis  for  a  computational  procedure  for  drawing  Monte  Carlo  sample  target 
paths.  The  mean  target  track  and  the  estimation  error  covariance  are 
approximated  as  sample  statistics  of  the  set  of  Monte  Carlo  tracks. 

As  we  have  posed  the  ADTR  problem,  there  is  a  time  (denoted  by  fj)  at 
which  the  initial  high  resolution  measurement  is  made  and  the  tracking 
problem  transitions  from  one  of  six-state  position /velocity  estimation  to  one 
of  full-state  (10-state)  estimation.  We  let  denote  the  time  of  the  last  low 
resolution  measurement  prior  to  f^.  We  refer  to  the  low  resolution  tracking 
phase  ending  at  tg  as  stage  1  and  to  the  joint  low  resolution  and  high 
resolution  tracking  phase  beginning  at  fj  as  stage  2. 

Standard  Kalman  filtering  methods  are  adequate  for  stage  1  tracking. 
The  assumption  then  is  that  at  time  tg  a  six-state  filter  estimate  of  the  target 
position  and  velocity  (and  hence  the  direction  of  target  flight)  is  available 
together  with  an  associated  6x6  error  covariance  matrix.  As  described  in 
Section  1.3,  the  stage  1  estimate  of  aircraft  direction  is  used  to  draw  statistical 
inference  about  the  target  aircraft  orientation  angles.  The  result  is  that  the 
stage  1  filter  output  provides  the  basis  to  construct  a  joint  distribution  on 
target  position,  velocity,  and  orientation  to  initiate  the  stage  2  (high 
resolution)  tracking  phase. 

The  state  estimation  technique  we  use  in  stage  2  is  fundamentally 
different  from  that  used  in  stage  1.  In  stage  2  we  use  a  form  of  Gibbs  sampling 
that  enables  random  draws  to  be  made  from  the  posterior  target  track 
likelihood  function  without  the  need  to  construct  the  full  likelihood  surface. 

The  particular  form  of  Gibbs  sampling  tihat  we  use  (see  reference  [3])  is 
based  on  the  Gibbs  representation  of  the  Bayes'  posterior  measure  Ptfdx): 


Pt(dx)  = 


(I-l) 
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The  function  H^  is  referred  to  as  the  potential  function  associated  with  |ij(dx). 
We  observe  that  up  to  an  additive  constant  H,  is  simply  the  negative 
log-likelihood  associated  with  the  Bayes'  posterior  measure. 

In  the  current  application,  is  defined  on  the  space  of  possible  target 
tracks  X(t)  on  [0,f]  sampled  at  a  sequence  of  discrete  times  .  In  general,  it  can 
then  be  shown  that  H,  has  an  additive  decomposition  of  the  form  = 
such  that  each  element  of  prior  information  about  X(f)  and  each  element  of 
information  about  X(f)  provided  by  a  particular  sensor  measurement 
separately  give  rise  to  a  summand  in  the  potential  function 
decomposition. 

The  component  likelihoods  which  contribute  to  H^iXit))  are 

•  the  prior  joint  likelihood  of  the  sequence  of  target 
position/ velocity  vectors  associated  with  the  track  X(f) 

•  the  prior  joint  likelihood  of  flie  sequence  of  target  Euler  angle 
vectors  associated  with  the  track  X(f) 

•  the  likelihood  of  measuring  the  target  position  at  time  to 
be  Zjt  given  target  track  X(f) 

•  the  likelihood  of  measuring  the  target  pixel  image  at  time  fjt 
to  be  fj  given  target  track  X(f). 


In  Section  1-3  we  will  derive  expressions  for  the  above  likelihood 
function  components  as  well  as  expressions  for  the  associated  potential 
function  components  p,.  In  the  case  of  position  and  velocity,  we  will 
postulate  a  stochastic  differential  equation  from  which  the  required  prior 
likelihood  can  be  derived.  We  will  also  postulate  a  Markov  model  for  the 
time  evolution  of  target  orientation.  In  this  case,  the  target  orientation  model 
will  define  one  likelihood  function  component.  A  second  component  will  be 
defined  by  the  constraint  that  the  target  pitch  and  yaw  angles  must  be 
consistent  (in  a  statistical  sense  we  will  make  precise)  with  the  direction  of  the 
target  velocity  vector. 

For  the  particular  sensors  modeled  in  our  Phase  I  study,  the  p,  arising 
from  sensor  measurements  can  all  be  dealt  with  within  the  common 
framework  of  measurements  with  Gaussian  error.  In  that  case,  the 
corresponding  potential  function  contributions  take  the  general  form 


p,-  (x)  =  (z  -  Zp  (x))^  (z  -  Zp  (x)). 
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where 


z  =  observed  measiirement  vector 
Zp(x)  =  predicted  measurement  vector  conditioned  on  target  state  x 
2'^  =  error  covariance  of  observed  measurement  vector. 


The  quantity  on  the  RHS  of  equation  (1-2)  is  the  so-called  Mahalonobis 
distance,  measuring  in  a  statistical  sense  the  distance  between  the  predicted 
measurement  Zp(x)  conditioned  on  the  assumed  target  state  x  at  the  time  of 
the  measurement  and  the  observed  measurement  z. 

In  the  case  of  a  low  resolution  sensor  measurement,  Zp(x)  is  the  target 
position.  Hence,  there  is  a  simple  analytic  relationship  (projection  onto  a 
subspace)  between  the  target  state  x  and  the  predicted  measurement.  The 
corresponding  potential  function  term  p,(x)  is  therefore  also  a  known 
function  of  x,  and  the  required  differentiations  in  calculating  the  gradient  of 
the  potential  function  can  be  carried  out  analytically. 

In  the  case  of  an  image  produced  by  the  high  resolution  optical  sensor, 
the  function  Zp(x)  depends  in  a  complex  way  on  the  target  state.  The  values 
of  Zp(x)  are  known  only  numerically  and  are  stored  in  a  precomputed 
template  library.  Consistent  with  our  model  for  the  optical  images  produced 
by  the  high  resolution  sensor,  each  template  image  is  represented  by  a  64  x64 
array  of  gray  scale  intensity  levels  (discretized  as  integer  values  on  the 
interval  [0,255]).  A  separate  template  is  stored  for  each  combination  of  aircraft 
orientation  and  aircraft  type.  For  this  purpose,  roll  angle  and  pitch  angle  are 
each  discretized  in  10  degree  increments.  The  step  size  in  yaw  angle  varies 
with  pitch  angle  in  such  a  way  that  the  density  of  grid  points  on  the  yaw-pitch 
sphere  is  approximately  uniform. 

All  of  the  stored  image  templates  are  based  on  a  reference  image  size 
scale  and  a  reference  azimuth-elevation  coordinate  origin  in  the  focal  plane 
of  the  imaging  sensor.  The  tracking  filter  generates  the  target  range  estimate 
to  rescale  each  observed  optical  image  to  the  reference  range.  We  assume  that 
background  pixels  can  be  distinguished  from  aircraft  body  frame  image  pixels 
and  take  the  center  of  area  of  the  projected  target  image  as  the  coordinate 
origin.  In  the  context  of  our  numerical  studies,  we  simulate  aircraft  pixel 
images  directly  in  the  template  library  viewing  frame  and  thereby  bypass  the 
image  registration  step. 

We  model  the  noise  in  each  pixel  as  additive  Gaussian  and  assume 
statistical  independence  from  pixel  to  pixel.  The  combination  of  the  template 
library,  the  Gaussian  noise  model,  and  simulated  image  data  provides  the 
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basis  for  determining  the  associated  with  each  high  resolution 
measurement. 

We  note  that  the  template  library  used  in  our  Phase  I  study  was 
generated  by  Dr.  Michael  Miller  and  his  colleagues  of  the  Department  of 
Electrical  Engineering  at  Washington  University  in  St.  Louis  (WUSL).  The 
images  in  the  library  are  all  for  a  standard  set  of  optical  viewing  conditions 
with  regard  to  atmospheric  clarity  and  background  light  intensity  level.  In 
Phase  n  we  plan  to  expand  the  template  library  to  include  a  range  of  viewing 
environments. 

The  jump-diffusion  method  proceeds  by  writing  down  a  stochastic 
differential  equation  (SDE)  of  the  Langevin  type  involving  the  gradient  VH, 
of  the  posterior  potential  fimction: 


dXit)  =  -^QQ^VHt(X(f))df  +  QdWit),  (1-3) 


where  dW(t)  is  a  standardized  white  noise  process  of  the  appropriate 
dimension  and  the  matrix  Q  governs  the  scale  of  the  diffusion.  In 
application  Q  is  taken  to  be  a  diagonal  matrix.  This  SDE  defines  a  continuous 
path  stochastic  process  describing  a  randomized  gradient  descent  on  the 
posterior  potential  function 

The  jump  element  of  the  model  comes  into  play  in  the  recognition  of 
the  target  type.  Target  type  is  a  discrete  variable,  and  a  mechanism  different 
from  the  continuous  track  deformation  implied  by  equation  (1-3)  is  required 
to  allow  transitions  in  target  type.  Specifically,  it  must  be  possible  for  the 
target  type  assigned  to  a  Monte  Carlo  sample  track  to  change  abruptly  from  a 
to  h  such  that  the  probability  of  such  a  change  is  determined  by  the  relative 
weight  of  statistical  evidence  favoring  b  as  the  target  type  as  opposed  to  a. 

The  particular  scheme  we  use  for  determining  when  such  discrete 
changes  in  target  type  occur  is  based  on  an  adaptation  of  the  Metropolis 
algorithm  (see  reference  [6]).  A  key  relevant  quantity  is 


H,(X(ft)la)  =  the  posterior  potential  of  X(t^)  given  target  is  of  type  a. 

The  decision  about  whether  to  alter  the  target  type  assigned  to  a  particular 
sample  track  is  based  on  H^iXit^^)Ha). 
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Suppose  that  the  target  type  assigned  to  the  sample  track  X(4)  at  time 
is  flg  and  that  the  alternative  target  type  possibilities  are  A  random 

draw  is  made  from  the  set  of  alternative  target  types  based  on  their  relative 
likelihood  of  occurrence.  Suppose  that  target  type  a^  is  drawn.  Then  the  rule 
for  assigning  the  new  target  t^e  b  to  X(t^)  is 


b  =  a,  if  H(X(f*)lfl,)  <  H(X(f,)lflo). 


(1-4) 


Otherwise, 


with  probability!- 


As  previously  noted,  it  is  shown  in  reference  [2]  that  the 
jump-diffusion  model  has  certain  ergodic  properties  which  insure  that  the 
process  has  a  limiting  stationary  distribution  equal  to  the  posterior 
distribution  on  the  target  track  conditioned  on  all  past  information  (both 
prior  information  and  sensor  data).  This  result  provides  the  theoretical  basis 
for  the  discrete  sampling  procedure  we  now  describe  for  estimating  the  target 
track. 


Let  f^,  k  =  l,2,”-,K,  denote  the  stage  2  measurement  times.  Then 
consider  the  stochastic  difference  equation  which  approximates  the  stochastic 
differential  equation  in  equation  (1-3): 


AX(f)  =  -^QQ'^VHt(X(t))Af +QAW(f).  (1-6) 


The  required  sample  paths  are  initialized  at  by  making  N 
independent  draws  from  the  joint  distribution  on  target  position  and 
orientation  based  on  the  time  tracking  solution  derived  from  the  output  of 
the  stage  1  tracking  filter.  The  assignment  of  a  target  type  to  each  path  is  also 
made  on  the  basis  of  Monte  Carlo  draws,  consistent  with  the  (assumed 
known)  frequency  fimction  for  occurrence  of  target  type.  We  note  that  our 
numerical  results  in  Chapter  n  are  based  on  N  =  100  sample  paths. 

Once  the  N  target  tracks  have  been  initialized,  the  next  operation  is  to 
extend  each  track  from  fp  to  fi*  Two  steps  are  involved:  (i)  the  forward 
projection  of  each  track  from  time  fp  to  time  fj  and  (ii)  the  updating  of  each 
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track  for  the  time  fj  measurement  data.  The  forward  projection  step  for  each 
path  is  based  on  a  dynamical  model  for  predicting  the  target  track  at  tj 
conditioned  on  the  sensor  data  through  This  dynamical  model  is 
described  in  detail  in  Section  1-3.1. 

The  measurement  update  at  time  is  accomplished  through  the 
application  of  the  stochastic  difference  equation  in  equation  (1-6).  An 
essential  element  of  the  theory  in  equation  (1-6)  is  that  it  must  be  iterated 
multiple  times  to  incorporate  the  effect  of  the  measurement  at  fj  into  the 
target  track  estimate.  Each  iteration  is  referred  to  as  a  cycle.  A  sufficient 
number  of  cycles  M  must  be  applied  so  that  the  target  track  arrived  at  under 
the  action  of  the  diffusion  equation  represents  a  random  draw  from  the  Bayes 
posterior  distribution  on  the  target  state  conditioned  on  all  of  the  past  data. 

The  same  two  step  process  is  applied  in  updating  the  sample  tracks 
from  time  f,-  to  time  First,  the  dead-reckoned  estimate  of  the  target  track 
at  is  determined  based  on  the  d)mamical  model  for  the  target  state.  Then 
the  jump-diffusion  model  is  applied  at  f,+|  to  incorporate  the  measurement  at 
^j+i  into  the  state  estimate. 

1-2.2  Computational  Considerations.  The  theory  provides  no  explicit 
guidance  on  how  large  M  must  be  at  each  f,  or  how  small  the  discrete  steps 
Af  at  each  cycle  must  be.  Our  experience  with  the  numerical  application  of 
the  diffusion  model,  however,  has  provided  the  following  insights: 

•  In  applying  the  diffusion  equation,  the  matrix  QQ^  in 
equation  (1-6)  effectively  scales  the  step  size  in  each  coordinate 
and  provides  the  option  to  choose  a  different  step  size  for  each  of 
the  three  categories  of  state  variable:  position,  velocity,  and 
orientation.  In  each  case.  At  is  chosen  so  that  the  change  in  state 
induced  by  a  single  diffusion  cycle  is  small  compared  to  the 
precision  with  which  the  mean  target  track  is  to  be  determined. 
The  numerical  results  in  Chapter  n  are  based  on  the  values: 


APmax  =  5  meters 

Ap„ax  =  1  meter  /  sec  (1-7) 

A<t)„,,  =  2.5°. 


In  each  case,  the  corresponding  At  is  obtained  by  substituting  the 
value  of  -^VH,(X(f))  for  the  current  cycle  into  equation  (1-6), 
assuming  the  contribution  of  the  white  noise  term  to  be 
and  solving  the  resulting  quadratic  in  At.  The  effect  then  is  that 
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At  is  dynamically  adjusted  for  each  category  of  state  variable 
following  each  cycle. 


There  is  a  tendency  for  a  very  large  number  of  cycles  to  be 
required  to  achieve  solution  convergence  in  equation  (1-6)  if  the 
posterior  mean  track  estimate  at  f,  differs  significantly  from  the 
prior  mean  track  at  f,-.  This  situation  occurs  when  a 
measurement  with  significant  corrective  power  (high 
information  content)  such  as  the  initial  high  resolution 
measurement  is  processed.  In  that  case,  it  becomes  necessary  to 
significantly  increase  the  tolerances  in  equation  (1-7)  to  accelerate 
convergence.  The  standard  tolerances  are  then  restored  for 
subsequent  update  times. 

In  obtaining  the  numerical  results  in  Chapter  II,  we  used 
M  =  2,000  cycles  except  in  processing  the  initial  high  resolution 
measurement  at  tj.  In  that  single  instance,  we  set  M  =  5,000. 


We  observe  that  the  principle  at  work  in  our  random  sampling 
technique  based  on  equation  (1-3)  involves  the  interaction  between  the 
gradient  of  the  potential  term  and  the  white  noise  term.  The  -^VH,  term 
tends  to  deform  the  sample  path  so  as  to  reduce  its  Gibbs  potential  and  hence 
to  increase  its  likelihood.  The  white  noise  term  injects  randomness.  This 
accomplishes  two  things.  First,  it  prevents  the  path  from  getting  stuck  at  a 
local  minimum  of  the  potential  function.  Second,  it  injects  the  correct 
amoxmt  of  randomness  into  the  diffusion  process  to  assure  that  the  cycling 
process  generates  sample  paths  approximately  in  proportion  to  their  relative 
likelihood. 

One  very  important  advantage  of  the  jump-diffusion  method  is  that 
the  sampling  of  the  likelihood  function  is  accomplished  without  having  to 
build  the  entire  likelihood  surface.  While  the  jump-diffusion  method  can 
require  significant  computation,  the  amount  of  computation  required  to 
construct  the  full  likelihood  surface,  because  of  the  high  dimensionality  of 
the  target  track  state  space,  would  be  prohibitive. 

The  application  of  equation  (1-3)  requires  the  calculation  of  the  gradient 
of  the  potential  function.  Because  of  the  additive  form  of  H,,  the  quantity 
-|VH,  can  be  expressed  in  terms  of  the  quantities  V  Analytic  expressions 
can  be  derived  for  all  of  the  V  p,  except  for  the  terms  arising  from  the  target 
pixel  images.  In  that  case,  the  are  specified  numerically  based  on  the  pixel 
images  in  the  template  library.  A  first-order  finite  difference  method  is  then 
used  to  calculate  each  V^,  . 


1-9 


In  Section  1-3  below,  we  will  give  explicit  expressions  for  each  of  the 
potential  terms  p,.  Let  x(/)  =  (p(/),t>(/))  denote  the  target  position/ velocity  at 
time  tj.  Then  one  term  of  the  potential  function  at  time  denoted 
^,(x(l),x(2),-  -,x(^C))  arises  from  the  joint  distribution  of  the  target  position 
and  velocity  at  the  observation  times  conditioned  on  the  past  data 

through  time  fx-i*  The  calculation  of  requires  knowledge  of  the 

covariance  matrix  for  x(K)  conditioned  on  the  data  through  The  sample 
covariance  for  x(fC)  obtained  from  the  M  Monte  Carlo  sample  paths  at  but 
prior  to  processing  the  observation  at  fj^  is  used  as  an  estimate  for  this 
covariance.  The  value  M  =  100  was  chosen  with  this  requirement  in  mind. 

We  observe  that  as  the  time  index  k  increases,  the  amount  of 
computation  required  to  execute  a  tracking  update  grows  at  an  approximately 
linear  rate.  One  way  of  inhibiting  this  growth  is  to  restrict  the  segment  of  the 
target  track  that  is  updated  to  a  moving  window.  If  the  length  of  this  window 
is  set  to  include  the  L  most  recent  high  resolution  measurements,  then  the 
measurement  at  time  would  be  used  to  modify  the  earlier  target  states  only 
at  those  earlier  measurement  times  tj  such  that  k- j <L-1.  The  result  of 
employing  such  a  moving  window  is  that  the  amount  of  computation  grows 
linearly  only  up  until  the  time  of  the  L*  high  resolution  measurement. 
Beyond  that  point,  the  computational  load  per  update  remains  roughly 
constant. 

The  results  in  Chapter  II  are  based  on  a  window  of  length  one.  Our 
estimation  methods  then  constitute  a  form  of  filtering  —  the  estimation  of 
the  current  system  state  conditioned  on  all  past  data.  We  note  that  the 
algorithm  code  has  been  written  to  handle  processing  windows  of  any  user 
specified  length.  In  its  more  general  mode  of  operation,  the  algorithm  can  be 
used  to  perform  the  smoothing  function  over  a  moving  window  of  fixed 
length  or  over  the  entire  track  history. 

We  observe  that  the  jump-diffusion  approach  to  sampling  from  the 
posterior  distribution  on  the  target  track  differs  fundamentally  from 
conventional  approaches  to  Monte  Carlo  tracking.  The  usual  procedime 
begins  by  postulating  a  prior  distribution  on  the  space  of  possible  tracks. 
Monte  Carlo  samples  are  drawn  from  this  prior.  These  sample  tracks  are  then 
assigned  relative  weights  based  on  the  likelihood  of  the  observed  sensor 
measurements  given  each  track.  As  more  and  more  measurements  are  used 
to  refine  the  sample  tracks,  the  most  unlikely  tracks  are  pruned  while  those 
that  survive  are  replicated  (the  regeneration  process)  to  insure  adequate 
statistical  coverage  in  the  shrinking  neighborhood  in  track  space  to  which 
ground  truth  becomes  progressively  confined. 
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Thus,  a  basic  feature  of  a  conventional  Monte  Carlo  tracker  is  that  the 
refinement  of  track  estimates  is  achieved  through  the  operations  of  pruning 
and  regeneration.  Individual  sample  tracks  are  extended  forward  in  time  but 
otherwise  do  not  change  in  response  to  new  data.  What  does  change  are  the 
relative  weights  assigned  to  the  sample  tracks. 

By  contrast,  in  the  jump-diffusion  model  approach  all  sample  tracks 
receive  the  same  relative  weight.  It  is  the  tracks  themselves  which  are 
deformed  in  response  to  the  data  through  a  combination  of  continuous 
variation  and  discrete  jumps.  It  is  this  deformation  process  which  governs 
the  statistical  properties  and  convergence  properties  of  the  sample  track 
estimates. 

The  tracking  and  identification  of  an  aircraft  leads  to  a  complex 
posterior  likelihood  function  and  a  high-dimensional  state  estimation 
problem.  In  such  an  application,  the  inherent  efficiency  of  the  path 
deformation  method  in  concentrating  the  Monte  Carlo  sample  paths  in  the 
relevant  region(s)  of  track  space  is  a  critical  factor  in  choosing  that  method 
over  the  path  re-weighting  method. 

We  now  turn  in  Section  1-3  to  a  more  detailed  specification  of  our 
ADTR  algorithm  than  we  have  attempted  to  tiiis  point. 


I>3  ADTR  Algorithm  Specification 

In  separate  subsections  below  we  describe  the  model  used  for  aircraft 
dynamics  and  the  stage  2  (high  resolution)  ADTR  algorithm. 


1-3.1  Aircraft  Dynamics.  We  begin  with  a  brief  description  of  our 
model  for  aircraft  dynamics.  Additional  details  are  provided  in  the 
Appendix. 

As  previously  described,  the  instantaneous  dynamical  state  of  an 
aircraft  is  represented  in  inertial  coordinates  in  terms  of  a  nine-vector 
comprised  of  a  three- vector  p-ip\,P2>V3^  specifying  the  aircraft  position,  a 
three- vector  u  =  (nj,n2,i>3)  specifying  the  aircraft  velocity,  and  a  three- vector 
<))  =  (<|)i,<j)2,<))3)  of  Euler  angles  specifying  the  aircraft  orientation  (<j)i  =  roll 
angle,  <1)2  =  pitch  angle,  <1)3  =  yaw  angle).  Roll  angle  and  yaw  angle  eacii  take 
values  on  the  interval  [0,27c]  with  the  endpoints  0  and  In  identified.  Pitch 
angle  takes  values  on  the  interval  [-n/2,n/7\.  The  aircraft  track  state  X(f)  at 
time  t  specifies  the  continuous  time  state  history  over  the  time  interval  [0,t] 
and  therefore 


X(f)  €  (91^  X  [0,27lf  X  [-71  /  2,71  /  2])'°''^. 
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The  particular  inertial  frame  we  choose  is  the  topocentric  coordinate 
frame  centered  at  the  common  location  of  the  low  resolution  sensor  and  the 
high  resolution  sensor.  The  positive  x-direction  is  local  east,  the  positive 
y-direction  is  local  north,  and  the  positive  z-direction  points  in  the  direction 
of  the  outward  normal  to  the  Earth  surface.  Azimuth  angle,  elevation  angle, 
and  range  are  the  spherical  coordinates  in  this  topocentric  frame.  In  the 
context  of  our  tracking  problem,  the  position  and  velocity  of  the  aircraft  are 
measured  relative  to  the  rotating  Earth.  This  allows  us  to  treat  the  Earth  as 
stationary  and  the  topocentric  frame  as  inertial. 

The  aircraft  orientation  angles  are  defined  by  the  sequence  of  axis 
rotations  defining  the  transformation  from  the  aircraft  body  frame,  which  has 
its  coordinate  axes  fixed  along  the  principal  axes  of  the  airframe,  back  into  the 
inertial  frame.  This  transformation  is  given  by 


‘1 

0 

0 

'cos(^2^ 

0 

-sm(<|)2)' 

0 

cos(<|)i) 

sinC^j) 
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0 
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sm(^2) 
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cos(<|>2) 
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cos(<))3) 
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sm(<j)3)  0 
cos((l)3)  0 
0  1 
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The  transformation  plays  an  important  role  in  the  use  of  the  image 
library.  The  templates  as  catalogued  are  for  a  standardized  viewing  geometry 
which  places  the  observer  at  (0,1,0)  (i.e.,  at  unit  distance  along  the  positive 
y-axis).  In  general,  if  an  aircraft  with  orientation  <j)  is  viewed  at  azimuth  0 
and  elevation  <p,  then  the  appropriate  orientation  <|)'  for  purposes  of  library 
look-up  satisfies  the  relation  =T^R(0,(p),  where  R(0,<p)  rotates  the  inertid 
frame  into  the  viewing  frame. 

1-3.2  Stage  2  -  High  Resolution  Tracking.  We  consider  the  stage  2  case 
in  which  there  is  a  current  time  to  through  which  the  processing  has 
proceeded.  We  let  t^,t2,"’>tj  denote  a  sequence  of  subsequent  times  and 
describe  the  procedure  for  extending  the  processing  to  tj.  We  note  that  not  all 
of  the  tj  need  be  measurement  times. 

We  describe  the  stage  2  procedure  for  the  general  case  where  the  joint 
distribution  on  target  position  and  velocity  at  time  tg  is  multivariate 
Gaussian.  At  the  start  of  stage  2,  tg  is  taken  to  be  the  time  of  the  last  low 
resolution  sensor  measurement  prior  to  stage  2.  High  resolution  tracking  is 
then  initialized  based  on  the  multivariate  Gaussian  track  solution  output  at 
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tg  by  the  stage  1  tracking  filter.  Throughout  the  subsequent  stage  2  processing, 
Monte  Carlo  methods  are  used.  The  role  of  tg  is  then  played  by  the  start  time 
of  the  time  window  over  which  the  sample  target  tracks  are  to  be  updated  for 
sensor  data. 

We  proceed  to  describe  the  calculation  of  the  various  terms 
contributing  to  the  potential  function.  We  consider  first  the  term  arising 
from  prior  information  on  the  aircraft  Euler  angles. 

The  approximate  values  (1>2  and  <j)3  of  the  initial  pitch  angle  and  initial 
yaw  angle  of  the  aircraft  can  be  inferred  from  the  direction  of  ^f^).  i^sume 
that  the  actual  initial  pitch  angle  is  constrained  to  the  interval  [^2  -e2/J>2 
and  the  actual  initial  yaw  angle  is  constrained  to  the  interval  [^3  -£3, <1)3 +  63]. 
No  information  is  assumed  known  about  the  initial  aircraft  roll  angle  other 
than  that  it  is  constrained  to  the  interval  [0,2jtl. 

Under  these  conditions,  the  probability  density  on  the  initial  triple  of 
Euler  angles  (<f>i/<|)2/^3)  is  taken  to  have  the  form  (for  n/2-l^2^»  82,83): 


C0S((|)2) 


87ce3Cos(<j)2  )sm(e2 )] 
0 


for  <|)j  e  [0,27t],  (j),-  e  [(J),.  - e, •,<!),  +  e,],  i  =  2,3 
otherwise. 


(1-9) 


We  use  the  Von-Mises  maximum  entropy  distribution  on  the  circle  to 
model  prior  uncertainty.  The  prior  conditional  density  on  <j)(/)  given  <})(/- 1) 
for  /  S  2  is  thus  assumed  to  have  the  Markov  Von-Mises  form 


1  K,cos[<|>.(y>».(;'-i)l  (T_|o) 

i|27tIo(Kf) 


Here,  k  =  (Ki,K2,K3)  is  a  vector  of  parameters  characterizing  the  degree  of 
concentration  of  Euler  angle  prior  density.  We  assume  that  each  k,.  is  of  the 
form  Kg(At)~^,  where  At  is  the  length  of  the  prediction  time  step.  The  scale 
parameters  Kq  for  roll,  pitch,  and  yaw  are  model  inputs. 
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It  follows  from  equation  (I- 10)  that  the  contribution  to  the  potential 
function  arising  from  the  Markov  transition  model  for  the  aircraft  Euler 
angles  has  the  form 


p,(<Kl),<l)(2),-,(l)(K))  =  -X  XK,cos((l),(y) (Ml) 


We  next  turn  to  the  potential  term  arising  from  the  position/ velocity 
model. 

We  assume  a  stochastic  differential  equation  model  for  aircraft  position 
and  velocity  of  the  particularly  simple  form 


dp  =  vdt 
dv  =  BdW(t), 


(1-12) 


where  dWit)  is  a  three-dimensional  white  noise  process  scaled  by  B  to 
produce  random  accelerations  with  covariance  matrix  BB^.  We  assume  that 
the  acceleration  process  is  isotropic,  so  that  BB^  has  the  special  form  o^I.  (In 
Phase  n  we  will  consider  the  case  of  non-isotropic  random  accelerations.) 

The  SDE  in  equation  (1-12)  is  linear  and  lends  itself  to  the  standard 
recursive  estimation  methods  of  Kalman  filtering.  Let  x(j)  =  ip{j),vij))  and 
assume  that  x(0)  is  given.  Under  the  action  of  the  dynamical  model  in 
equation  (1-12),  the  conditional  distribution  for  x(;  +  l)  given  x(;)  is 
multivariate  Gaussian.  For  /  =  0,1,  •vfC-l  define 


li;+i  =  E[x(/  +  l)lx(y)] 

^;+i  =  Cov[x(  ;■  + 1),  x(; ■  + 1)1  x( ;)]. 


Then  py+j  =  Ayx(/),  where  the  matrix  Ay  can  be  expressed  in  terms  of  the 
coefficients  in  equation  (1-12).  The  Kalman  filter  prediction  equations  also 
provide  a  recursive  relationship  for  covariance  matrices  Ey. 

It  follows  that  the  potential  function  contribution  associated  with 
the  joint  distribution  of  the  x(/),  ;  =1,2,  ",K,  can  be  decomposed  as 
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(1-14) 


p,(x(l),x(2),-,x(K))=i(x(7)-|i,fsf(x(/)-4,). 


It  is  an  easy  matter  to  calculate  Vp,  from  the  above  expression. 

The  next  potential  term  arises  from  the  assumed  relationship  between 
the  direction  of  the  aircraft  velocity  vector  and  the  aircraft  orientation  angles. 
Let 


0(/)  =  azimuth  angle  of  vij) 


=  tari 


-1 


(p(/)  =  elevation  angle  of  p(;) 


=  tari 


-1 


(1-15) 


(Note:  We  assume  that  the  aircraft  velocity  vector  provides  no  information 
about  the  aircraft  roll  angle.) 

We  then  let 


Pi  =  -Xf^el  COs((pij)  -  <|>2(/))  +  COs(Q(j)  -  (1)3(/))],  (1-16) 

/=! 


where  and  are  constants  chosen  to  fit  the  Von-Mises  density  to  the 
density  g(<l)i,<t>2/<l>3)  equation  (1-8).  We  note  that  the  effect  of  this  p,  term 
on  the  overall  potential  function  is  mathematically  equivalent  to  that  of  a 
measurement  of  the  direction  of  the  aircraft  velocity  vector  given  the  aircraft 
pitch  and  yaw  or,  equivalently,  a  measurement  of  the  aircraft  pitch  and  yaw 
given  the  direction  of  the  aircraft  velocity  vector. 

Monte  Carlo  sample  tracks  are  extended  from  to  in  two  stages. 
The  first  stage  is  to  randomly  draw  the  orientation  angles  at  4+^.  This  is  done 
using  the  Von  Mises  density  in  equation  (I-IO)  governing  orientation  angle 
transitions.  The  second  stage  is  to  randomly  draw  the  aircraft  position  and 


1-15 


velocity  at  fj^+i  based  on  equations  (1-12)  and  (1-13).  Once  the  Monte  Carlo 
tracks  have  been  projected  forward  to  the  diffusion  model  is  applied  to 
incorporate  the  effect  of  sensor  data. 

A  summary  of  the  steps  involved  in  the  execution  of  stage  2  of  the 
ADTR  algorithm  are  as  follows: 


Step  1.  Initialize  N  Monte  Carlo  tracks  at  time  to  by  making  random 
draws,  first  from  the  aircraft  position/velocity  density  (based 
on  the  low  resolution  tracking  filter)  and  then  from  the 
conditional  density  on  the  aircraft  orientation  (equation  (1-8)). 
Set  jt  =  0. 

Step  2.  Project  the  time  Monte  Carlo  tracks  forward  in  time  to 

based  on  the  SDE  model  for  position/velocity  and  the 
Von-Mises  model  for  orientation.  Set  k  =  k  +  l. 

Step  3.  Process  the  sensor  measurements  at  times  using 

equation  (1-3)  and  a  random  sampling  method  based  on  the 
jump-diffusion  model.  Execute  M  sampling  cycles  for  each 
Monte  Carlo  track  X(4),  reevaluating  VHj(X(f|t))  after  each 
cycle.  Return  to  step  2. 


The  above  procedure  proceeds  iteratively  until  all  of  the  scheduled 
updates  are  executed  or  until  there  are  no  more  measurements  to  be 
processed.  After  the  measurement  update  is  completed  at  each  f^,  the 
following  sample  statistics  are  calculated  and  output: 


(i)  the  sample  means  for  the  aircraft  position  vector,  velocity 
vector,  and  orientation  vector 

(ii)  the  three  submatrices  of  the  error  covariance  matrix 
corresponding  to  position,  velocity,  and  orientation 

(iii)  the  fraction  of  Monte  Carlo  paths  assigned  each  of  the  possible 
target  t5q)es. 


In  our  Phase  I  implementation  of  the  algorithm,  we  assume  the  target 
aircraft  to  be  of  known  type  and  focus  on  the  tracking  function.  The 
numerical  results  we  present  in  Chapter  II  compare  the  sample  mean 
estimate  for  target  position  and  the  (simulated)  ground  truth  target  position. 
The  magnitude  of  the  vector  difference  between  the  estimated  and  ground 
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truth  target  positions  is  the  assumed  measure  of  tracking  error.  Our  interest 
then  is  in  how  tracking  error  behaves  as  a  function  of  time,  particularly 
during  a  target  maneuver. 

We  note  that  the  algorithm  we  have  implemented  in  Phase  I  differs 
from  that  described  in  reference  [3]  in  certain  important  respects. 

The  algorithm  described  in  reference  [3]  tracks  the  target  in  the  six-state 
space  of  position  and  orientation.  We  have  added  velocity  to  the  target  state 
vector.  Since  the  translational  motion  of  the  aircraft  is  a  Markov  process 
jointly  in  position  and  velocity,  the  state  prediction  step  can  then  be 
performed  recursively  in  much  the  same  manner  as  in  Kalman  filtering. 
Furthermore,  when  velocity  is  explicitly  part  of  the  target  state,  the 
relationship  between  the  target  (translational)  direction  of  motion  and  the 
target  orientation  can  be  enforced  directly  through  the  Gibbs  potential 
function. 

The  algorithm  in  [3]  merges  initial  detection  and  the  initiation  of  high 
resolution  tracking  into  a  single  integrated  operation.  The  point  of  view  we 
have  taken  is  that  because  the  high  resolution  sensor  has  a  narrow  field  of 
view,  in  an  operational  setting  it  is  likely  that  track  initiation  would  be 
performed  by  a  low  resolution  search  sensor  with  a  much  wider  field  of  view. 
High  resolution  tracking  would  then  be  initiated  via  a  direct  handoff  from 
the  low  resolution  sensor.  We  have  taken  this  two-stage  approach  in 
modeling  the  tracking  process.  A  computational  benefit  is  that  it  eliminates 
the  computation  intensive  function  of  having  to  search  through  the  entire 
template  library  to  initiate  high  resolution  tracking. 
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Chapter  II 


Test  Studies 


In  this  chapter  we  summarize  the  results  of  our  numerical  studies 
applying  the  tracking  algorithm.  In  Section  11- 1  we  define  the  test  cases  we 
analyzed.  In  Section  11-2  we  present  numerical  results  and  in  Section  11-3 
discuss  our  Phase  I  study  conclusions. 


II-l  Test  Case  Descriptions 

The  test  cases  we  ran  were  designed  with  two  objectives  in  mind.  First, 
we  wanted  to  quantify  the  benefits  of  multisehsor  tracking.  In  particular,  we 
wanted  to  demonstrate  that  the  error  in  tracking  a  target  aircraft  can  be 
reduced  significantly  by  fusing  high  resolution  optical  image  data  with  low 
resolution  radar  data.  Second,  we  wanted  to  investigate  the  robustness  of  the 
tracking  accuracy  of  the  algorithm  to  variations  in  the  characteristics  of  the 
high  resolution  sensor  as  well  as  to  variations  in  the  resolution  and  accuracy 
of  the  image  template  library. 

We  first  describe  the  assumed  ground  truth  aircraft  track.  We  then 
identify  the  assumptions  made  about  the  sensors  and  the  image  library, 

II-1.1  Target  Aircraft  Track.  In  general,  we  define  a  ground  truth 
aircraft  track  in  terms  of  the  following  parameters: 

•  initial  aircraft  position  vector  p(0)  and  orientation  angle 
vector  <j>(0) 

•  aircraft  speed  v  (assumed  to  be  constant) 
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•  a  sequence  of  time  intervals  i  =  (withfo=0) 

and  orientation  angle  rate  vectors  <j>(i)  such  that  ^{i)  is  the 
assumed  constant  aircraft  orientation  angle  rate  vector 
on  [f.vfoil. 


We  assume  that  at  any  instant  the  direction  of  aircraft  flight  is  the 
positive  x-direction  in  the  body  coordinate  frame.  As  a  consequence  of  this 
assumption,  the  above  data  are  sufficient  to  determine  the  full 
nine-dimensional  position,  velocity,  and  orientation  state  of  the  aircraft  as  a 
function  of  time. 

The  assumed  aircraft  track  as  shown  in  Figure  II-l  consists  of  an  initial 
level  flight  leg,  a  banked  and  pitched  turn  of  180°,  and  a  second  level  flight 
leg  as  shown  in  Figure  II-l.  The  aircraft  is  initially  at  10,000  ft  and  ascends  at  a 
constant  pitch  of  20°  while  in  the  turn.  The  aircraft  roll  angle  in  the  turn 
is  -45°.  The  aircraft  yaw  rate  (°/sec)  during  the  turn  is  denoted  ^3.  We  note 
that  ^3!?  is  the  magnitude  of  the  acceleration  implied  by  the  aircraft  yaw 
rate  ^3  and  speed  v.  We  have  taken  ^3  =  5°/sec,  so  that  the  180°  turn  requires 
36  seconds  to  complete.  A  500  mi  /  hr  aircraft  speed  and  a  5° /sec  yaw  rate  in 
the  turn  corresponds  to  approximately  a  2g  aircraft  acceleration.  The  aircraft 
altitude  upon  completing  the  turn  is  approximately  19,000  ft. 

II-1.2  Sensor  Characteristics.  The  characteristics  of  the  low  resolution 
sensor  are  representative  of  those  of  a  military  air  traffic  control  radar  or  a 
shipboard  radar  and  are  treated  as  fixed^.  The  characteristics  of  the  high 
resolution  optical  sensor  are  varied  parametrically.  The  principal  value 
shown  below  for  each  parameter  represents  the  base  case.  Alternate 
variations  are  shown  in  parentheses. 


^  Such  radars  scan  in  azinauth  and  do  not  measure  target  elevation  angle.  However,  for 
present  purposes  we  have  given  the  low  resolution  radar  a  minimal  capability  (Oq  =  2°)  to 
measure  elevation  angle. 
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Ground  Truth  Aircraft  Track 


IT-3 


Low  Resolution  Sensor 


Measurement  errors: 

Oaz  =  0.5« 

Gq  =  2.0° 

Or  =  200  meters 


Tracking  rate: 


Plr  =  0.2  Hz 


High  Resolution  Sensor 

Measurement  errors: 

aA,  =  0.r  (0.05°,0.2°) 
00  =  0.1°  (0.05°,0.2°) 


Tracking  rate^: 


Phr  =  2Hz  (IHz,  4  Hz). 


II-1.3  Image  Templates.  The  image  template  library  for  the  X29  was 
provided  by  WUSL.  The  template  resolution  is  64  x  64  pixels  with  the 
reflected  visible  light  intensity  level  represented  on  a  gray  scale  indexed  from 
0  to  255.  Roll  and  pitch  angles  are  discretized  imiformly  in  10°  increments. 
The  yaw  angle  increment  varies  with  pitch  angle.  The  number  of  yaw  values 
is  maximum  at  32  (a  step  sizell.25°)  when  the  pitch  angle  is  0°.  The  yaw 
discretization  becomes  increasingly  coarse  as  pitch  increases  until  there  is 
only  a  single  yaw  value  of  0°  at  a  pitch  of  90°. 

Observed  images  are  simulated  by  adding  (truncated)  Gaussian  white 
noise  to  the  library  templates.  The  scale  of  the  image  noise  is  measured  in 
terms  of  a  number  of  gray  scale  intensity  levels.  The  assumed  base  case  pixel 
noise  level  is  o^  =  12  gray  scale  units  (gsu). 


2  The  actual  frame  rate  of  the  optical  sensor  may  be  20  Hz  -  60  Hz.  However,  we  assume  a 
much  lower  effective  tracking  rate  to  reflect  any  multi-frame  integration  inherent  in  the 
signal  processing  and  also  to  justify  our  assumption  that  angle  measurement  errors  are 
statistically  independent  and  that  the  pixel  noise  in  the  optical  images  is  white. 
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The  base  case  parameter  values  and  variations  used  in  the  image 
simulation  and  processing  are  summarized  as  follows: 

pixel  resolution  =  64  x  64  (32  x  32,  16  x  16) 
gray  scale  levels  =  256  (128,  64) 

pixel  noise  =  12  gsu  (24  gsu,  48  gsu). 


As  shown  in  Figure  II-l,  the  two  sensors  are  collocated  at  a  point  which 
is  50  miles  (measured  horizontally)  from  the  aircraft  at  the  midpoint  of  its 
turn.  Low  resolution  tracking  begins  at  f  =  0.  The  high  resolution  sensor 
begins  reporting  at  f  =  20  secs,  and  from  that  point  on  tracking  is  based  on 
fusing  the  low  resolution  and  high  resolution  data.  The  aircraft  180°  turn  is 
initiated  at  t  =  30  secs  and  completed  at  f  =  66  secs,  at  which  point  the  tracking 
problem  is  terminated. 

We  analyzed  three  tracking  variations  based  on  which  sensor  data  are 
fused  into  the  track  solution: 


(LR);  Low  resolution  position  measurements  only 

(HRA):  Low  resolution  position  measurements  +  high 

resolution  angle  measurements 

(HRI):  Low  resolution  position  measurements  +  high 

resolution  angle  measurements  +  high  resolution 
pixel  images. 


The  track  estimation  in  the  LR  and  HRA  processing  variations  is  based 
on  the  Kalman  filter  described  in  Chapter  I.  In  the  HRI  variation,  the  image 
data  is  added  and  the  jump-diffusion  model  approach  (also  described 
in  Chapter  I)  is  used. 


II-2  Numerical  Results 

Figures  II-2  through  II-ll  below  represent  the  principal  numerical 
results  of  our  Phase  I  study.  Figures  n-2  through  II-6  apply  to  the  base  case 
and  show  in  succession:  error  in  position,  velocity,  roll  angle,  pitch  angle,  and 
yaw  angle.  In  the  case  of  position  error  and  velocity  error,  the  LR  and  HRA 
processing  results  are  shown  for  comparison. 
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The  position  tracking  results  in  Figure  11-2  show  that  during  the  period 
from  t  =  20  secs  to  f  =  30  secs,  while  the  aircraft  is  still  in  level  flight,  the 
fusion  of  the  high  resolution  angle  data  with  the  low  resolution  position  data 
has  a  very  significant  effect.  The  tracking  error  rises  to  as  high  as  3,000  meters 
if  only  the  radar  data  are  processed.  If  the  optical  sensor  angle  measurements 
are  fused  into  the  tracking  solution,  the  tracking  error  is  kept  under 
500  meters. 

Once  the  aircraft  maneuver  begins  at  f  =  30  secs,  even  the  combination 
of  radar  data  and  optical  angle  data  is  not  enough  to  keep  the  tracking  error 
from  increasing  from  about  200  meters  at  the  beginning  of  the  maneuver  to 
as  much  as  1,400  meters  during  the  maneuver.  One  observes  that  in  this 
HRA  processing  case  the  position  error  exhibits  a  very  distinct  sawtooth 
behavior  with  a  5  sec  period.  The  position  error  drops  sharply  each  time  the 
radar  reports  a  target  range  measurement  and  rises  systematically  between 
successive  radar  measurements.  This  behavior  reflects  the  fact  that  the 
optical  sensor  provides  no  direct  target  range  information  and  therefore  the 
component  of  position  error  along  the  line  of  sight  tends  to  grow  between 
radar  measurements. 

Once  the  image  data  are  fused  into  the  tracking  solutions,  the  situation 
changes  rather  markedly.  Now  the  target  orientation  information  deduced 
from  the  target  images  can  be  used  to  deduce  the  direction  of  the  target 
forward  motion.  As  a  result,  knowledge  of  the  target  position  along  the  line 
of  site  does  not  degrade  between  radar  observations.  One  observes  in 
Figure  II-2  that  the  sawtooth  behavior  prominent  in  the  HRA  case  is  now 
much  less  in  evidence.  The  net  result  is  that  position  error  under  HRI 
processing  is  consistently  less  than  under  HRA  processing  and  remains  in  the 
range  200  meters  to  600  meters,  about  the  same  as  during  the  period  the  target 
was  in  level  flight. 

Figure  n-3  shows  the  base  case  velocity  (vector)  error  comparison. 
Again  the  benefit  of  using  imagery  data  in  support  of  tracking  is  in  evidence. 
During  the  interval  from  f  =  20  secs,  to  f  =  30  secs.,  while  the  target  is  still  in 
level  flight,  velocity  error  is  consistently  less  than  30  m/sec  under  HRA 
processing.  Once  the  target  starts  its  maneuver,  there  is  a  more  or  less 
systematic  increase  in  velocity  error  until  about  40  secs,  both  under  HRA  and 
HRI  processing.  This  early  stage  of  the  maneuver  presents  a  particularly 
difficult  tracking  problem  since  two  and  sometimes  all  three  of  the  target 
orientation  angles  are  changing  simultaneously.  Once  the  target  roll  and 
pitch  steady  to  their  new  values  and  only  the  target  yaw  angle  is  changing, 
velocity  error  trends  generally  downward  under  HRI  processing.  Under 
HRA  processing,  velocity  error  trends  generally  upward  imtil  it  peaks  at  about 
200  m/sec  at  f  =  50  secs,  (the  approximate  midpoint  time  of  the  target 
maneuver)  before  starting  a  slow  decline. 
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Figures  II-4  through  II-6  show  the  behavior  in  the  base  case  of  the 
errors  in  estimates  of  the  target  orientation  angles  imder  HRI  processing.  One 
observes  that  yaw  angle  error  tends  to  be  larger  in  magnitude  and  more 
variable  than  either  roll  angle  error  or  pitch  angle  error.  One  contributing 
factor  is  that  once  the  target  maneuver  starts,  the  target  yaw  angle  varies 
continuously  at  the  rate  of  5° /sec  throughout  the  36  sec.  period  of  the 
maneuver.  Roll  and  pitch,  on  the  other  hand,  stabilize  at  a  constant  value 
within  9  secs.  There  is  also  the  effect  that  the  discretization  of  yaw  angle  is 
variable  with  pitch  angle  and  generally  coarser  than  the  discretization  of  roll 
and  pitch. 


We  summarize  our  base  case  numerical  results  as  follows: 


Case  I: 

(i)  target  not  maneuvering 

(ii)  fusion  of  low  resolution  (radar)  and  high  resolution 
(optical)  point  tracking  data  (HRA  processing) 

(iii)  Kalman  filter  state  estimation  methods. 

Then  tracking  errors  vary  from  200  meters  to  500  meters.  This 
compares  with  tracking  errors  of  as  much  as  3,000  meters  under 
low  resolution  (LR)  point  tracking. 


•  Case  II: 

(i)  target  maneuvering  (2g  acceleration  for  36  secs.) 

(ii)  fusion  of  low  resolution  point  tracking  data,  high 
resolution  point  tracking  data,  and  imaging  data  (HRI 
processing) 

(iii)  jump-diffusion  method  of  state  estimation. 

Then  tracking  errors  vary  from  200  meters  to  600  meters.  This 
compares  with  tracking  errors  of  as  much  as  1,400  meters  under 
high  resolution  (HR)  point  tracking. 


Base  Case:  Roll  Error 
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Figures  TL-7  through  II-ll  show  the  effect  on  tracking  error  of  varying 
five  of  the  key  problem  parameters  relative  to  their  base  case  values: 

•  optical  sensor  direction  measurement  error 

•  optical  sensor  measurement  rate 

•  number  of  pixels  per  target  image 

•  number  of  pixel  intensity  levels 

•  pixel  noise  level. 


In  each  case,  tracking  error  is  plotted  out  to  the  midpoint  of  the  aircraft 
maneuver  (f  =  48  secs).  Each  figure  shows  three  plots:  the  base  case  and  two 
variations  of  the  indicated  parameter. 

Figures  11-7  and  n-8  show  that  there  is  a  marked  sensitivity  of  tracking 
error  to  the  data  rate  and  measurement  error  of  the  optical  sensor.  One 
observes,  in  particular,  that  a  factor  of  two  variation  either  in  data  rate  or 
direction  measurement  error  can  induce  a  factor  of  two  variation  in  tracking 
error. 


Figure  11-9  shows  the  effect  of  varying  the  number  of  pixels  used  to 
represent  an  aircraft  image.  We  note  that  as  the  number  of  pixels  is 
decreased,  an  adjustment  is  made  in  the  pixel  noise  level  to  accovmt  for  the 
aggregation  of  smaller  pixels  into  larger  pixels.  For  example,  each  pixel  in  a 
32  x  32  image  is  the  composite  of  4  pixels  in  the  corresponding  64x64 
representation,  such  that  the  intensity  level  assigned  to  the  composite  pixel  is 
the  average  of  the  intensities  of  the  4  component  pixels.  This  averaging  has 
the  statistical  effect  of  reducing  the  pixel  noise  level  by  a  factor  of  two. 

One  observes  from  Figure  n-9  that  the  relationship  between  the  level 
of  pixel  aggregation  and  tracking  error  appears  to  be  a  complex  one.  Prior  to 
about  36  secs,  tracking  error  is  minimal  when  32x32  pixel  images  (as 
opposed  to  16  X 16  pixel  images  or  64  x  64  pixel  images)  are  used.  Beyond 
about  42  secs,  tracking  error  is  minimal  when  64  x  64  images  are  used. 
Between  36  secs  and  42  secs,  the  three  pixel  aggregation  levels  result  in  about 
the  same  tracking  error. 

A  possible  explanation  for  this  observed  behavior  is  that  the  optimal 
degree  of  pixel  aggregation  depends  on  the  target  viewing  geometry.  More 
pixel  aggregation  implies  less  image  detail  but  reduced  noise  in  each  pixel. 
Less  aggregation  implies  more  image  detail  but  increased  noise  in  each  pixel. 
The  best  operating  point  is  likely  to  vary  with  the  amount  of  structural  detail 
of  the  target  that  is  needed  to  resolve  its  orientation.  This  in  turn  varies  with 
the  target  viewing  geometry. 
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Measurement  Error  Variations 
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Data  Rate  Variations 
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Variations  in  Number  of  Pixels 
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Figure  n-10  shows  the  result  of  decreasing  the  number  of  distinct  gray 
scale  levels.  As  in  the  pixel  aggregation  case,  it  is  appropriate  to  adjust  the 
pixel  noise  level  as  the  niraiber  of  gray  scale  levels  is  varied.  In  the  present 
case,  the  pixel  error  was  adjusted  to  maintain  a  constant  ratio  between  the 
noise  level  measured  in  gray  scale  units  and  the  number  of  gray  scale  levels. 

Figure  II-IO  shows  that  tracking  error  consistently  increases  as  the 
number  of  gray  scale  intensity  levels  decreases.  Furthermore,  the  effect  is 
typically  a  rather  pronounced  one. 

Figure  II-ll  shows  the  effect  of  varying  the  level  of  pixel  noise.  One 
observes  that  there  is  something  of  a  threshold  effect  at  work,  in  that  the 
increase  in  noise  from  12  gsu  to  24  gsu  has  only  a  minimal  effect  on  tracking 
error.  However,  the  further  increase  in  pixel  noise  from  24  gsu  to  48  gsu  has  a 
marked  impact  in  terms  of  increased  tracking  error. 

We  summarize  our  observations  based  on  our  parametric  study  as 
follows: 

•  Tracking  error  shows  a  strong  dependence  on  the  measurement 
accmacy  and  data  rate  of  the  optic^  sensor 

•  A  decrease  in  the  number  of  gray  scale  intensity  levels  and  an 
increase  in  the  level  of  pixel  noise  each  cause  the  tracking  error 
to  increase.  The  magnitude  of  the  effect  is  small,  however,  in 
the  case  of  increasing  the  pixel  noise  from  12  gsu  to  24  gsu 

•  The  effect  of  varying  the  number  of  pixels  per  image  appears  to 
be  a  complex  one,  suggesting  a  possible  geometry-dependent 
tradeoff  between  image  detail  and  effective  pixel  noise  level. 


We  have  made  no  systematic  effort  in  Phase  I  to  benchmark  the  fusion 
algorithm  execution  time.  Optimization  of  the  algorithm  implementation  in 
order  to  meet  operational  execution  time  constraints  will  be  a  Phase  II 
technical  objective.  We  do  observe,  however,  that  the  execution  time  per 
measurement  update  varies  approximately  linearly  with 


N  =  (number  of  Monte  Carlo  sample  tracks) 

X  (number  of  diffusion  cycles  per  measurement  update). 
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Variations  in  Number  of  intensity  Leveis 
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Pixel  Noise  Variations 
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In  our  Phase  I  testing,  we  used  N  =  100x2,000  =  200,000,  a  value 
deliberately  chosen  generously  high  from  the  standpoint  of  the  statistical 
convergence  of  track  estimates.  In  Phase  II  we  will  develop  improved 
methods  for  choosing  the  step  size  per  cycle  in  an  effort  to  reduce  significandy 
the  number  of  diffusion  cycles  required.  We  will  then  investigate  in  detail 
the  trade-offs  between  the  statistic^  accuracy  of  the  algorithm  estimates  and 
execution  time.  The  eventual  operational  goal  in  Phase  HI  will  be  to  achieve 
near  real-time  processing  on  the  designated  computing  platform. 


II-3  Study  Conclusions 

Our  principal  Phase  I  study  conclusions  are: 

•  We  successfully  demonstrated  that  the  computations  required  to 
carry  out  the  data  fusion  algorithm  operations  using  our 
estimation  methods,  previously  executed  only  within  the 
context  of  a  massively  parallel  computer  architecture,  can  be 
performed  on  a  serial  processing  workstation  (the  equivalent  of 
a  Sim  Microsystems  SPARC  20). 

•  Our  parametric  study  shows  that  there  are  significant  benefits  to 
be  gained  in  terms  of  increased  tracking  accuracy  through  the 
fusion  of  multisensor  point  tracking  and  target  image  data.  We 
have  demonstrated  this  specifically  for  the  combination  of  a  low 
resolution  radar  sensor  and  a  high  resolution  optical  imaging 
sensor.  However,  the  methods  we  have  developed  to  fuse 
multisensor  data  are  quite  general  and  can  be  extended  to  other 
sensing  technologies. 
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Appendix 


Aircraft  Dynamics 


In  this  Appendix  we  record  certain  formulas  related  to  the  d5mamics  of 
aircraft  motion.  These  formulas  are  useful  in  constructing  ground  truth 
aircraft  tracks  and  in  determining  the  accelerations  implied  by  hypothesized 
aircraft  maneuvers. 

Consider  the  aircraft  body  frame  which  has  the  coordinate  axes  fixed 
along  the  principal  axes  of  the  airframe  and  define 


=  aircraft  velocity  vector  projected  onto  body  frame 

=  aircraft  angular  velocity  vector  projected  onto  body  frame 
/b  = 

=  vector  of  forces  acting  on  aircraft  in  the  body  frame 


Rigid  body  d)mamics  requires  that  the  following  system  of  differential 
equations  must  hold  at  any  time  s  (see  Section  2.6  of  [7]): 


■  cil(s)vl{s)  +  qlis)vlis)=  /f(s) 
62(5)+  ^3(s)z?f(s)  -  qf(s)??3(s)=  /f(s) 
UgCs)  -  ^/2(s)^^f(s)  +  <lf(s)P2(s)  =  /“(s). 


(A-1) 


We  note  that  (A-1)  can  be  expressed  in  vector  form  as  /b  =  +  ^b/ 

where  the  term  x  represents  the  Coriolis  acceleration  induced  by  the 
rotation  of  the  body  frame  relative  to  a  reference  inertial  frame. 

The  rotational  velocity  vector  ^^7  frame  is  equal  to 

the  vector  (<i>i,<j>2/^3)/  where 

=  roll  rate 
^2  =  pitch  rate 
^3  =  yaw  rate. 

The  angular  velocity  (9i/<72/%)  f*'  ffr®  inertial  frame  is  related  to  the  Euler 
angles  by 

<7i(s)  =  <j)i(s)  -  <i>3(s)sm((t)2(s)) 

^2<s)  =  <i>2(s)cos(<t)i(s)) + ^3(s)cos((|>2(s))sm(<t),(s))  (A-2) 

q^is)  -  -^2fs)sin(<j)i(s))+<j>3(s)cos(<>2(s))cos(<|>i(s)). 

The  inertial  position  vector  is  given  by 

pis)  =  r  ^'^ix)v^ix)dx  +  pito),  (A-3) 

where  the  transformation 

'10  o'  "cos(<l)2(x))  0  -si«((|>2(x)) 

T(t)=  0  cos(<|)i(x))  sm(<l)i(t))  X  0  1  0 

0  -sini^iix))  cosi^-^ix))  sini^2^x))  0  cos((|)2(x)) 

’  cosi^^ix))  sini^3ix))  0 
X  -sm(<j>3(x))  cos(<l)3(x))  0 
0  0  1. 

rotates  the  inertial  frame  axes  into  the  body  frame  axes. 


A-2 


f' 


If  one  assumes  that  the  aircraft  velocity  vector  in  the  body  frame  points 
in  the  positive  x-axis  direction,  then  =  {v,0,0),  where  v  is  the  aircraft  speed. 
If  p  is  assumed  constant,  then  the  relationship  /b  =  ^  can  be  used  to 

relate  a  specified  combination  of  roll  rate,  pitch  rate,  and  yaw  rate  to  the 
implied  accelerations.  It  then  becomes  a  straightforward  matter  to  design 
aircraft  maneuvers  involving  accelerations  of  a  specified  number  of  g’s, 
where  g  is  tiie  gravitational  acceleration.  If  the  assumed  Euler  angle  rates  are 
assumed  constant,  one  can  directly  perform  the  integration  in  Equation  (A-3) 
to  determine  the  aircraft  position  as  a  function  of  time. 
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